#ifndef VEC3_H_
#define VEC3_H_

#include "../common.h"

#include "Vec.h"

namespace physicsworks {

namespace math {

using exception::EUndefMathOpt;

template<class T> class Vec3 : public Vec<T> {
private:
	typedef Vec<T> Base;

public:
	T &x, &y, &z;
	
	using Vec<T>::operator*;
	using Vec<T>::operator();
	using Vec<T>::operator-;

	Vec3<T>() : Vec<T>(3), x(Base::arr[0]), y(Base::arr[1]), z(Base::arr[2]) {}
	Vec3<T>(const Vec<T>& v) : Vec<T>(v), x(Base::arr[0]), y(Base::arr[1]), z(Base::arr[2]) {
		// The provided vector must be a three dimensional one
	}
	Vec3<T>(const Vec3<T>& v) : Vec<T>(v), x(Base::arr[0]), y(Base::arr[1]), z(Base::arr[2]) {}
	Vec3<T>(T X, T Y, T Z) : Vec<T>(3), x(Base::arr[0]), y(Base::arr[1]), z(Base::arr[2]) {
		x = X; y = Y; z = Z;
	}
	Vec3<T>(const T v[3]) : Vec<T>(3), x(Base::arr[0]), y(Base::arr[1]), z(Base::arr[2]) {
		x = v[0]; y = v[1]; z = v[2];
	}

	/**
	 * Returns the vector's length.
	 */
	T length() const {
		return sqrt((*this) * (*this));
	}

	/**
	 * Normalizes the vector.
	 * 
	 * To deal with very small vectors, it first scales the vector by 1/abs(max(x,y,z))
	 * then normalizes it.
	 */
	 Vec3<T>& normalize() {
	 	std::valarray<T> arra = std::abs(Base::arr);
	 	*this /= arra.max();
	 	T l = length();
	 	
	 	Assert<EUndefMathOpt>(!DMATH || !equal(l, T(0)), "The vector's length must be nonzero to normalize it");
	 	
	 	return *this = (*this) / l;
	 }

	/**
	 * Normalizes the vector. Old version.
	 */
//	void normalize() {
//		T l = length();
//		
//		Assert<EUndefMathOpt>(!DMATH || !equal(l, T(0)), "The vector's length must be nonzero to normalize it");
//		
//		*this = (*this) / l;
//	}

	/**
	 * Returns the dot product.
	 */
	const T operator * (const Vec3<T>& v) const {
		T r = 0;
		for (size_t i = 0; i < 3; i++) {
			r += this->arr[i] * v.arr[i];
		}
		return r;
	}
	
	const Vec3<T> operator - (const Vec3<T>& v) const {
		Vec3<T> r = *this;
		return r -= v;
	}
	
	/**
	 * Assigns the other vector.
	 */
	Vec3<T>& operator = (const Vec3<T>& v) {
		Base::operator=(v);
		return *this;
	}
	
	/**
	 * Returns the cross product.
	 */
	Vec3<T> cross(const Vec3<T>& v) const {
		Vec3<T> r;

		r.x =  y * v.z - z * v.y;
		r.y = -x * v.z + z * v.x;
		r.z =  x * v.y - y * v.x;

		return r;
	}
	
	void operator () (const Vec<T>& v) {
		x = v(0);
		y = v(1);
		z = v(2);
	}
	
	void operator () (T X, T Y, T Z) {
		x = X; y = Y; z = Z;
	}
	
	static Vec3<T> X;
	static Vec3<T> Y;
	static Vec3<T> Z;
	static Vec3<T> ZERO;
};

template<class T> Vec3<T> Vec3<T>::X    = Vec3<T>(1, 0, 0);
template<class T> Vec3<T> Vec3<T>::Y    = Vec3<T>(0, 1, 0);
template<class T> Vec3<T> Vec3<T>::Z    = Vec3<T>(0, 0, 1);
template<class T> Vec3<T> Vec3<T>::ZERO = Vec3<T>(0, 0, 0);

}

}

#endif /*VEC3_H_*/
